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Expressions are derived for the indefinite integrals, 

jrf(r) C (ar)dr 

jrj(r) Q(ar)dr 

(rf(r) C (ar) Co(fir)dr a^j3 

J //(/) Co(ar) Zo(kr)dr 

where Co{ar) are zero order Bessel functions, Zo(kr) are zero order modified Bessel functions and /V) is a polynomial in r.- 
In general, the expressions given for the integrals are given in terms of prescribed functions of the Bessel functions, and the 
coefficients of these functions are determined from a finite series, the terms of which are found from recurrence relation- 
ships that involve only the polynomial /(r). Coefficients of the terms of the finite series are given in tabular form for up to 
an eleventh degree polynomial. 
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1. Introduction 

In this paper four indefinite integrals are studied; namely, 

\rf(r) C (ar)dr (1) 

jrf(r) C%(ar)dr (2) 

jrf(r) Co(ar) C ((3r)dr a^($ (3) 

jrf(r) Co(ar) Z (kr)dr, (4) 

where f(r) is a p\h degree polynomial in r, 

f(r) = J a n r\ 

Co(ar) = AJ {) {ar) + BY {) {ar) , and Z {) (kr) = DI {) (\r) + EK (\r) . Jo(ar) and Yo(ar) are zero order 
Bessel functions of the first and second kind, respectively, and Io(kr) and Ko{kr) are zero order 
modified Bessel functions of the first and second kind, respectively. Expressions for integrals 
(1), (2), and (3) can be found in Watson [l] 1 and Luke [2] which must be solved individually for 
every term of the polynomial /(r) by rather cumbersome recursion formulas involving the Bessel 
functions. The expressions given here for the integrals are given in terms of prescribed functions 



1 Figures in brackets indicate the literature references at the end of this paper. 
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of the Bessel functions, and the coefficients of these functions are determined from a finite series, 
the terms of which are found from recurrence relationships that involve only the polynomial /(r). 
Butler and Pohlhausen [3] have provided solutions of (1), (2), and (3) as definite integrals that 
specifically apply to solutions where Jo(a) =Jo(/3) = 0. A more general application is needed. 
A method, not hitherto available in the literature, is presented here to enable a user to determine 
these indefinite integrals involving Bessel functions. Further comment in evaluation of the method 
is given at the end of the paper. 

An example of the need for solutions to the integrals (1-4) is illustrated in the study of physical 
problems involving the steady-state diffusion process in a solid of cylindrical geometry, where a 
particular partial differential equation to be satisfied is 

d 2 v 1 dv , d 2 v 
c)r- r dr dz 2 

v is the potential, and r and z are the radial and longitudinal dimensions, respectively. A solution 
of (5) for the cylinder of finite length and outer radius is of the form 

fsinh az] „ fsin kz] 

v = C {) (ar)\ \+Zo(\r)\ . }, (6) 

[cosh az) [cos kz) ' 

where, for the solid cylinder B = E = 0, or C {) (ar) = AJo(ar) and Z () {kr) = D/o(Xr), and final solution 
of (6) is dependent upon the evaluation of A. B, Z), and E from the boundary conditions on the 
curved and plane surfaces of the cylinder. For rather simple boundary conditions, the constants 
can readily be determined by use of the orthogonality properties of the Bessel and trigonometric 
functions. 

More difficult boundary conditions on one or more of the plane and curved surfaces of the 
cylinder may take the general forms 

civ 

—=H(r) [v—V(r)\ on a plane surface (7) 

az 

civ 

— =h(z) [v—V{z)\ on a curved surface (8) 

dr 

where the proportionality factors H(r) and /?(z), and the potentials of the ambient V(r) and V(z) are 
position dependent. To obtain values for the constants of (6) upon the substitution of (6) and its 
first derivatives in (7) and (8), it would be convenient to use relationships such as 

k cos X + h(0) sin A = (9) 

and aCt(a)-H(0)C (a) = (10) 

and their related orthogonality properties, and positive roots. For the determination of the constants 
A and B of (6), expressions are necessary for the integrals (1-4) with appropriate limits of integration 
assigned. 

The reason for the specific boundary conditions cited above is to further and to ameliorate 
the study of radiant heat transfer and other nonlinear boundary conditions over nonisopotential 
surfaces in conjunction with heat conduction and other diffusion processes. 



of frf(r) 



2. Evaluation of \rf(r) Co(ar)dr 

From expressions for the integral (1) for the special cases /(r) =r m , m= 1, 2, 3, and 4, derived 
by integration by parts, it can be shown by induction that a general expression for (1) is of the form 

frf(r) C (ar)dr= L |"f(r) C, (or) + F ' (r) Ct>{ar) ] -E1M \c () {ar)dr (la) 
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where the polynomial F(r) defined by this relation is to be determined in the subsequent analysis, 
and 

&(ar)=AJ l (ar)+BY 1 (ar), 
\C (ar)dr= r€ (ar) +| [C,(ar) H (ar)-Co(ar) H,(ar)]r (lb) 



and //o and H, are Struve functions of the first kind. Differentiation of (la) gives 

F „ (r) + F , (r)-F , (0) +a2F(r) = Q2/ , (r) _ 



(11) 



The differential equation (11) is satisfied by the following relationships: 

n=2 



where for n=\ 



and for n > 1 



P l (r)=f"(r) + 



/'(r)-/'(0) 



P«(r) =P H n . l (r) + 



/>;,.,(>•) -^,,(0) 



In terms of the coefficients «,„ of a pth degree polynomial 

.S'=l 

Table 1 shows values of /{ m , n for values of m from 1 to 11. For the sake of example, let/(r) = r 7 , 
then 



p/ v 7 49r* 1225^ ll,025r 

F(r)=r 7 — — —-i — — -H: — , 



a L a* 



and 



r // x ™ 245/^ , 3675^ 11,025 
F (r) = 7H> — I — 

rul ^4 ^wb 



are determined for substitution in (la). 

Table 1. Values of R m>n 



V 












-\ 
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2 


3 


4 


5 


1 













2 


4 











3 


9 











4 


16 


64 









5 


25 


225 









6 


36 


576 


2,304 







7 


49 


1225 


11,025 







8 


64 


2304 


36,864 


147,456 





9 


81 


3969 


99,225 


893,025 





10 


100 


6400 


230,400 


3,686,400 


14,745,600 


11 


121 


9801 


480,249 


12,006,225 


108,056,025 
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263-048 O - 67 - 6 



An interesting sidelight is the substitution of Z (Xr) for Co(ar) in (la), which gives 
jrf(r)Z n (kr)dr= ~ [AZ,',(Xr)F, (r) -Z (Xr)F;(r)] +^^ jz (\r)dr 
where 

F.(r)-/(r) + 2^ 

and P w (r) is as defined for (la). 

3. Evaluat 



Hon of I rf(r)Cl(ar)dr 



/,/< 



An expression for (2) can be written in the form 
r)C%{ar)dr= r ^\G(r)^ 



+ 



2a 2 



G"(r) + -G'(r) 
r 



CUar) 



4a* J C,) 



(12) 



(ar)^r (2a) 



where the polynomial G(r) is to be determined in the subsequent analysis. The integral I J 2 (ar)dr 
is given in (23), for k=l. Differentiation of both sides of (2a) gives 



rG'{r) + G{r)=f{r) 



1 
4a 2 



rG"'(r) + 3G"(r)" 



G»-G'(0) 



The expression on the left can be written rG'{r) + G(r) = — [rG(r)\ then integration gives 



CW -IJ /W *--LJ 



rG"'(r) + 3G"(r) + 



G'(r)-G'(0) 



dr. 



(13) 



(14) 



G(r) is satisfied by 



rM v (~l)"g»(r) 

G < r »=S — ^ — 

n = " 



(15) 



where for n = 



^ 



So(r)=- \f(r)dr 



and for n > 



S " (r) = ir/ rS »-i( r ) + 3S »->( r > + 



S^Jr)-S^M 



dr. 



Also, (14) is satisfied for G(r) in terms of the coefficients a,,, of a pth degree polynomial by 



^here 



1 A (^-25-h2) 3 
r m>0 =l/(m+l),andr m , tt =^ TT n ( w -25+l) 



(15a) 
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Values of T mj „ are given in table 2 for m from to 1 1. For sake of example, let /(r) =ao + air 
-f a>r 2 + a;i7 3 + <7 4 r 4 , then from (15a) and table 2, 

8<7 2 . 27(/;sr . 64a 4 r 2 
^/i/ . (i 2 r 2 , </;{/-'* . ^/4/ 4 3 8 IS 512a 4 



G(r)=oo+^+-^+-|-- 



+ - 



(2«)* 2 15(2c*) 



where the first and seeond differentials of G(r) are necessary for substitution in (2a). 

Table 2. Values o/T m n 



v= 














\ 





1 


2 


3 


4 


5 


" \ 
















1 


1 
1/2 












2 
3 


1/3 

1/4 


8/3 

27/8 










4 


1/3 


64/15 


512/15 








5 


1/6 


125/24 


1425/16 








6 


1/7 


216/35 


1,608/35 


36.864/35 






7 


L/8 


343/48 


42.875/192 


385.875/128 






H 


1/9 


512/63 


12.288/35 


262.144/35 


2,097.152/35 




9 


1/10 


729/80 


83.349/160 


2.083.725A28 


56.260.575/256 




10 


1/11 


1000/99 


512.000/693 


2.457.600/77 


52.428.800/77 


419,430,400/77 


II 


1/12 


1331/120 


323.133/320 


37.979.173/640 


924.479.325/5 1 2 


24.960,941,775/1024 



4. Evaluation of rf{r)Co(ar)C {) ((3r)dr, a^ fi 
An expression for (3) may take the form 

[ rf(r)Co{ar)Co{(ir)dr=N x U i (r)+N 2 U\(r)-y\N;U-Ar) + N 4 U: i (r) \ 
- [i/»t/',(0) ~ y 2 Ui(0)] J C (ar)Co(j3r) dr + y[U[ (0) - <jtU! z (0) | J C , <«r)C , (/3r),/r 



t^here 



8 = a 2 -j8 2 , </> = (a 2 + £ 2 )/8 2 , y=2oc(ilb\ 



(3a) 



£/i(r) and f/a(r) are polynomials to be determined in the subsequent analysis, and 



N, 



- j rC„(ar)C l) (8r)dr = ^[aC l (ar)C„(l3r)- /3C n (ar)C,(l3r)] 



N, = j rC , (ar)C , (/3r) dr = ^ \(5C , (ar)C tt (fir) - aC tt (ar)C , (/3r) ] 

/V 2 /r = ^/Co(ar)C (/3r) + y C , (ar)C , (j8r) 
W 4 /r = 4iCt(ar)C , (/3r) + y C»(ar)C (j3r) 

The integrals I J t) (ar),l»((3r)dr and I J t {ar) J i((3r) dr are found in section 6. Differentiation of both 
sides of (3a) gives the following equations: 
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i/,(r) + i// 

U 2 (r) + ip 



UM+ 



U'M-U'M 



U&)- 



im-u'M 



UV) Ufr)-U'M 



-fif) 



= 0. 



Letting 



and 



A" = R = 



(16) 



A- = w = 

then substitution in (16) gives 
L ,o(r)=f(r) 



M k ,n(r) = L" k>n (r) 



LlnirtlLkjM 



+ M' 



k, n-l 



M' k . n -M)-M' k . n _M 

r 



(17) 



where, negative subscripts denote zero quantities. In terms of the coefficients a m of a pth orde 
polynomial, the quantities of (17) have been evaluated for m from zero to 11 in table 3. 



Table 3 









Coefficients of a m 


for y° 


Coefficients of a„, for y 2 


Coefficients 


of a,„ for y 4 


m 


r m 


r m-l 

Mo,o 


Mo,, 


r m 6 
Lo,3 

/Wo. 2 


U 4 
Mo, 3 


r m-\0 
Mo, A 


/-i.o 


M,, 


r m 8 

U , 
M ul 


Mi.i 


L 2 .0 


r m-10 
L>. 1 

M 2 ,„ 



1 

2 


1 



1 



1 




4 























3 


1 




9 
3. 






















4 


1 




16 
8 


64 









32 













5 


1 




25 

15 


225 
120 








135 













6 


1 




36 
24 


576 
480 


2,304 







384 



3,456 











7 


1 




49 

35 


1,225 
1,260 


11,025 

7,445 






875 



19,215 

2,625 










8 


1 




64 
48 


2,304 
2,688 


36,864 
39,936 


147,456 





1,728 



70,656 
13,824 


442,368 





55,296 





9 


1 




81 

63 


3,969 
5,040 


99,225 
135,135 


893,025 
703,080 




3,087 



203,175 
46,305 


3,044,790 
748,440 




416,745 





10 


1 




100 
80 


6,400 
8,640 


230,000 
360,960 


3,686,400 
4,730,880 


14,745,600 



5,120 



495,360 
122,880 


13,701,120 
4,945,920 


73,728,000 



1,966,080 



27,648,00 




11 


1 




121 
99 


9,801 
13,860 


480,249 
828,135 


12,006.225 
19,625,760 


108,056,025 
94,895,955 


8,019 



1,072,071 
280,665 


47,505,150 
20,291,040 


604,178,190 
203,388,570 


7,016,625 



245,768,985 
21,049,875 
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^(^/W-^oj + ^.u-. . . + y-[/,,,„-M///.,,, + t// 2 L,,2 • • .] + y 4 [^2,o-i//L-i//Z> 2 , 1 

+ . . .1 (18) 
f/2(r) = Mo,o-i//M ,, + i// 2 M () ,2- . . • +y 2 [Afi,o-<M#i,i+ • • .] + y 4 [W2,o-<M/2,i+ . . .]. (19) 

From the table we get 

f(r)=r* C/ 1 (r)=r2-4i/i 
£/2(r) = 

/(r)=r 3 £/,(r) = r 3 -9n// 
U 2 (r) = 3r 

f(r) = r 4 £/,(r) = r 4 -16r 2 i// + 64i// 2 + 32y 2 
£/ 2 (r) = 8r 2 



f(r)=r* f/ 1 (r) = r>-25r*i// + 225n// 2 + 135ry 2 
£/ 2 (r)=15r 3 -120n// 



/(r) = r« t/i(r)= r 6 - 36/ 4 i// + 576r 2 i// 2 - 2304i// 4 + 384y 2 r 2 - 3456y 2 i/; 
U 2 (r) = 24r 4 -4«0r 2 i/i 

/(r) = a () + air-h a 2 r 2 -f a.^r 3 U\(r) = a {) + a x r+ a 2 r 2 + oaf 8 — (4a 2 + 9a 3 r)i// 

t/ 2 (r)=3a 3 r. 



-/ 



5. Evaluation of J rf(r)Cn(ar)Zu(kr)dr 

An expression for (4) may take a form similar to (3a) 
j rf(r)C u (ar)Z u (kr)dr=W i V l (r)-W i V' 1 (r)+y[fV 3 V 2 (r)-W 4 V'(r)] 

+ ^V[(0) + wW: 2 (0)] |c (ar)Z„(Xr)^r- W [F;(0)+^(0)] f C 1 (ar)Z l (\r)dr (4a 

where c/> = a 2 + X 2 , f- (X 2 -a 2 )/</> 2 , <o = 2aX/4> 2 

r 1 =frCo(ar)Z (Xr)dr=^[aC,(ar)Zo(Xr) + XCo(ar)Z 1 (Xr)] 

r 3 =frC,(ar)Zi(Xr)dr=^[XCi(ar)Z (Xr)-aCo(ar)Zi(Xr)] 

— =fCo(ar)Zo(Xr)+a)C,(ar)Z,(Xr) 
r 

JF 4 

— = |Ci ( ar ) Zi ( Xr ) — coCo ( ar) Z ( Xr) 
r 

and Zi(\r)=DIi(kr)-EKi(kr). 
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An analysis similar to that of section (4) gives 

K,(r) = 2 (-l)*a**2fZ*.,(r) 

A:=0 n=0 

F 2 (r) = 2 (-l)*o) 2fc X^ Mfc .»( r ) 

k=Q n=0 

where Lk, n {r) and Mk,n(r) are as defined in section 4 and table 3. 

6. Evaluation of I Mar)M/3r)dr and I Ji{ar)J \{&r)dr 

The expression for the first integral may be written in the form 

I Jo(ar)J (l3r)dr= j J (ar)dr-a* j rf(r)J (ar)dr 

where from (la), the integral becomes 

f Mar)MPr)dr=ll+F'(0)] j Jo(ar)dr- r[aF(r)Mar)F' (r)J (ar)] 

and where 



f(r)=k? 



l -MM 

j8 2 r 



(r /8V gV 

\2 2 2 2 -4 2 2 2 -4 2 



*-•) 



(20) 

(21) 



is an infinite series and k = (5ja< 1. From section 2, F(r) becomes 



- (-l)"^'"- /- 2 "-' " (i) 2 m^ 2m 

(r) ^ ri 3 5 • • (2rc-l)l 2 ^ 



{2n-\)Y± {m\Y 



where (i),„ = r(i+m)/r(i). 

The complete elliptic integral of the first kind is defined 



K{k2)= L [l-A- 2 sin 2 ^]>/ 2= 2A 



The summation on /n in the above expression for F(r) is then defined by 

? h-1 (1^2 1.2m 



7T 



(m!) 2 



The integral upon substitution becomes 
j ./o(ar)./„(/3r)^/r = | /^(/,- 2 ) f 7 (ar)rfr 



\ n^v2(n-l) r 2n— 1 



+ 1 [1-3-5. , {2n _ V)] 2 ^(/r){ar./ l (ar) + (2n-l)J l) (ar)} 



(22) 



138 



where the integral f J {) (ar)dr has heen defined in section 2. The series is convergent for k< 1. 
An expression for the integral / Jo(fir)Y {) {(xr)dr is found by substituting Y {) {ar) , Y x (ar) for Jo(ar) 
and ,/i(ar), respectively, in (22). An expression for the integral / J { )(ar)I {) {/3r)dr is found from 
(22) by substitution /- 2 = - /r. //-'" = (- 1 X''"'. 

The expression (22) is limited for numerical evaluation due to a large loss of significance for 
large arguments although for values of/, 2 <^ 1 there is probably no restriction for numerical evalua- 
tion. Another more general expression for the integral can be obtained using the method of Butler 
and Pohlhausen [3, pp. 11-18], An integral representation of the Bessel function is 

>s (krt)dt 2 { n, ~ cos (r cos 0) sin Odd 



. ,, , 2 P cos (krt)dt 2 f^ 2 cos (r 

77 Jo (l-^ 2 ) 1 / 2 7T Jcos-1 A" (/<~ 



cos 2 0) 



1/2 



where the expression on the right is derived from the substitution, At = cos 0. A definite integral 
upon appropriate substitution becomes 

[ V Jo(ar)J (Mdr=- \" Jo(r)J Q (kr)dr 

J() OL J,) 

= — —T7 y „ -TT7Z ./()(/") COS (/COS 9) dl'(W 

no) cosi k (fr — cos- 6/) 1 /- Jo 
= y[7 2 U)+2./ 2 (z)]+2y./,(z) f t-'V^.U) y 



+ y./„(z) f (-DVtod) [^ 



L+ JW 
I 2n + I 



(231 



where z = ay, and 



_ = 2_ p/ 2 sin (2/i+ l)ft</6> 
2 " + ' 7rJ cos ,, (/r-eos 2 ^)'^' 



From [3, eq 55] values of ^2» + i are determined by the recurrence formula 

> = 2//,(2/r-l) _^ll^ , 

m + 1 /?? -f 1 

where ¥i = l, ^ 3 = 2A: 2 - 1, * 5 =6A: 4 -6A 2 + 1, etc. The integral for J 2 Q {ar)dr is found from (23) by 
setting k= 1, whereby S&2n+\ = 1. 

Numerical values for the definite integral 



/: 



Jo(ar)J (Pr)dr 



have been determined from (22) and (23) and are shown plotted against k for various values of a 
in figure 1. For a few cases, the numerical values have also been obtained by numerical integration 
by Simpson's rule. 

The expression for the second integral may be written by integration by parts in the form 

| J,(ar)J i (Mdr=- Mar) ^ r) +k j Mar)Mfir)dr-k f J -^-Jo(ar)dr. 

139 



1.0 




I(a,k) = J Jo(ar)Jo(/3r)dr 
Jo 






1 ' 1 ' 1 ' 1 
a=l 
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.8 










______2__ 


- 


.6 






-^^- 


Ka.iO 




— — — — ^— — — — _ 


- 


4 






-=: 




___a=4_^ 


— 


.2 




_J^^^^ ____--— 


i 




1 


10 — — -- "" 


n 


1.1.1,1 



.2 4 .6 .8 1.0 

k = /3/a 
Figure 1 

The first integral on the right has been solved previously ((22) or (23)), and the second integral from 
(la) takes the form 



(MM , 



Pr 



| + F'(0) I Mar)dr-r[aF(r)Mar)+F'(r)Mar)] 



nth fir) =P M2 - MM =P 



/3 2 r 3 | ffV 5 



2 2 • 4 2 2 • 4 2 • 6 2 2 • 4 2 • 6 2 



From section 2, F(r) becomes 



^ (-l)n a 2(n-l) r 2n-l , (j)jk^ . 

1 2, [1 • 3 • 5 • • ' (2/1- l)] 2 £„ (™!) 2 (2m + 2) 



An integral of K(k 2 ) is defined by 



f 



k'K(k' 2 )dk' = E(k'' 



(1 \2 L2m+2 



2,^ () ( TO !) 2 (2m + 2) 



where E(k 2 ) is the complete elliptic integral of the second kind. The summation on m in the above 
expression for F(r) is then defined by 



9 n-\ ( A\2 1.2m 

o„ ( F) =- [*(*) - a- W w] - g o ( J; 2 ;^ +2) 
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The integral upon substitution becomes 
f Ji(fSr)Ji(ar)dr=- Mar) J l{fir) +^ [*(*») ~E(k 2 )] j Jo(ar)dr 

+ / ' X n ( 7 1 ) w ^ w "' )r2w " 1 /t /l (/-2){ ar y l(ar ) + (2 B -1 )J Q (ar)} (24) 
iT, U '3 -5 • • ■ (2n-l)\- 

where 

n-l (O m 4 1W.U2 i.2m 

TTfr £ Q (2m+2)(m!) 2 

7. Discussion 

An example of the usefulness of the method presented in this paper is given for the evaluation 
of the definite integral 



I! 



rf(r)J (ar)dr (25) 



where /(r) = ao H~ fl^r 2 + a^r 3 -+- a.*/* 4 . The integral is of the form given in (la), where the coefficients 
of the polynominal F(r) are given in table 1. The resulting expressions for the upper limit are 

a 2 a 4 

f (1)=/ (1) , and 

a- 

F'(0)— 2a- 

From (la), the definite integral (25) becomes 

rf{r)Jo(ar)dr= 1 ; 1 - Jo{ar)dr, 

Jo ' a a- (X 4 J„ 

where the integral on the right is evaluated using (lb). Given values for a , a>>, a 3 , a 4 , and a, and 
tables of Bessel and Struve functions, the value of the definite integral is readily computed by desk 
calculator. If subroutines are available for computing the Bessel and Struve functions, digital 
computer calculation is rapidly accomplished. 

In this instance, the value of the definite integral is determined with need only for evaluation 
of F(r), using table 1, and then of F(l), F'(l), F'(0), and of four additional functions, namely, 
7o(«), Ji(ot) , Ho(a) and H\ (a). The computational task is minor as compared with that of evalu- 
ating many subdivisions or steps of the original integrand for a numerical integration of equivalent 
accuracy. 
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